# -*- coding: utf-8 -*-
"""
Created on Sat Apr  1 21:30:44 2017

@author: GangTimes
"""
import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D 
import matplotlib.font_manager as fm
#回归系数 [a0-an]
#因变量 [x0-xn] x0=1
#自变量 y={0,1}
#数据规模 
num=200
#绘图设置
zh_font=fm.FontProperties(fname=r"c:\windows\fonts\simsun.ttf",size=18)#路径用小写  大写不行 有这个不能保存为eps
en_font=fm.FontProperties(family='monospace',size=18)
zf_font=fm.FontProperties(family='Times New Roman',style='italic',size=16)

fig=plt.figure(figsize=(7,5))
ax=Axes3D(fig,[0.,0.0,1,1])
plt.rcParams['axes.unicode_minus']=False #用来正常显示负
ax.set_xlabel('x',fontproperties=zf_font)
ax.set_ylabel('y',fontproperties=zf_font)
ax.set_zlabel('z or p',fontproperties=zf_font)

def map_fun(theta,X):
    '''
    Logistic回归模型的极大似然函数
    极大似然函数是一个凹函数
    凹函数的极大值问题是凸优化问题
    此处采用了极大似然函数进行取负
    
    theta:回归系数 matrix[mx1]
    X:自变量（包含常数项）  matrix[nxm] n样本数 m代表自变量的维数
    Y:因变量({0,1})  matrix[nx1]
    return:matrix[1x1]
    '''
    zz=numpy.matrix(numpy.ones((num,1)))
    XX=numpy.hstack((X[:,0],X[:,1],zz))

    var_true=(theta.T)*XX.T*X[:,2]
    var_all=(numpy.log(1+numpy.power(numpy.e,(theta.T)*XX.T)))*zz

    return var_all-var_true

def create_data():
    '''
    用于生成一个二维空间的两种类型的点
    人为制造一定的误差
    return:matrix[nx3] n为样本数量  [x y z]三列
    '''
    x=numpy.random.uniform(0,1,num)
    y=numpy.random.uniform(0,1,num)
    z=list()
    for index in range(num):
        temp=numpy.random.random()
        xtemp=x[index]
        ytemp=y[index]
        if(temp>0.05):
            if(xtemp+ytemp>1):
                z.append(1)
            else:
                z.append(0)
        else:
            if(xtemp+ytemp>1):
                z.append(0)
            else:
                z.append(1)            
    
    x=numpy.matrix(x).T
    y=numpy.matrix(y).T
    z=numpy.matrix(z).T
    self_plot3d(x,y,z)
    return numpy.hstack((x,y,z))
    

def jacobi(theta,X):
    '''
    构建雅克比向量
    X:自变量构成的矩阵 matrix[nxm] n为样本数 m为参数个数[x,y,1]
    theta:参数构成的矩阵 matrix[mx1]
    return:matrix[mx1]
    '''
    x=X[:,0]
    y=X[:,1]
    z=X[:,2]
    zz=numpy.matrix(numpy.ones((num,1)))
    XX=numpy.hstack((X[:,0],X[:,1],zz))
    val_true_x=numpy.dot(x.T,z)
    val_true_y=numpy.dot(y.T,z)
    val_true_z=numpy.dot(zz.T,z)
    
    temp=numpy.power(numpy.e,XX*theta)/(1+numpy.power(numpy.e,XX*theta))
    
    val_all_x=numpy.dot(temp.T,x)
    val_all_y=numpy.dot(temp.T,y)
    val_all_z=numpy.dot(temp.T,zz) 

    return numpy.vstack((val_all_x-val_true_x,val_all_y-val_true_y,val_all_z-val_true_z))

def hessian(theta,X):
    '''
    构建hessian矩阵
    X:自变量构成的矩阵 matrix[nxm] n为样本数 m为参数个数[x,y,1]
    Y:样本的输出值
    theta:参数构成的矩阵 matrix[mx1]
    return:matrix[mx1]
    '''
    
    x=X[:,0]
    y=X[:,1]
    zz=numpy.matrix(numpy.ones((num,1)))
    XX=numpy.hstack((X[:,0],X[:,1],zz))
    temp=numpy.power(numpy.e,XX*theta)/(1+numpy.power(numpy.e,XX*theta))
    a11=numpy.sum(numpy.multiply(temp-numpy.power(temp,2),numpy.power(x,2)))
    a12=numpy.sum(numpy.multiply(temp-numpy.power(temp,2),numpy.multiply(x,y)))
    a13=numpy.sum(numpy.multiply(temp-numpy.power(temp,2),x))
    a21=a12
    a22=numpy.sum(numpy.multiply(temp-numpy.power(temp,2),numpy.power(y,2)))
    a23=numpy.sum(numpy.multiply(temp-numpy.power(temp,2),y))
    a31=a13
    a32=a23
    a33=numpy.sum(temp-numpy.power(temp,2))
    '''
    下面的写法不对
    a11=numpy.dot(temp.T,numpy.power(x,2))-numpy.dot(numpy.power(temp,2).T,numpy.power(x,2))
    a12=numpy.dot(temp.T,numpy.multiply(x,y))-numpy.dot(numpy.power(temp,2).T,numpy.multiply(x,y))
    a13=numpy.dot(temp.T,x)-numpy.dot(numpy.power(temp,2).T,x)
    a21=a12
    a22=numpy.dot(temp.T,numpy.power(y,2))-numpy.dot(numpy.power(temp,2).T,numpy.power(y,2))
    a23=numpy.dot(temp.T,y)-numpy.dot(numpy.power(temp,2).T,y)
    a31=a13
    a32=a23
    a33=numpy.dot(temp.T,z)-numpy.dot(numpy.power(temp,2).T,z)
    '''
    return numpy.matrix(numpy.vstack((numpy.hstack((a11,a12,a13)),numpy.hstack((a21,a22,a23)),numpy.hstack((a31,a32,a33)))))


def step_length(theta,delta,X):
    '''
    采用回溯直线搜索求取符合wolfe条件的步长
    theta:待估计参数 matrix[mx1]
    delta:newton步径 matrix[mx1]
    X:自变量构成的矩阵 matrix[nxm] n为样本数 m为参数个数[x,y,1]
    Y:样本的输出值 matrix[nx1]
    t:步长matrix[1x1]
    return: t 优化后的步长matrix[1x1]
    '''
    alpha=0.2
    beta=0.8
    t=1
    while(True):
        fvall=map_fun(theta+t*delta,X)
        fvalr=map_fun(theta,X)+alpha*t*jacobi(theta,X).T*delta     
        if(fvall[0,0]>fvalr[0,0]):
            t=beta*t
        else:
            break
    return t

def newton_main():
    '''
    newton 主算法程序
    '''
    X=create_data()
    error=0.001
    theta=numpy.matrix([[10,-0.5,20]]).T
    for i in range(100):    
        g=jacobi(theta,X)
        h=hessian(theta,X)
        delta=-(h.I)*g        
        decrease=(g.T)*(h.I)*g
        print(decrease)
        if(abs(decrease[0,0])<error):
            break
        t=step_length(theta,delta,X)
        theta=theta+t*delta      
    print(theta)
    logistic(theta)
    
def logistic(theta):
    '''
    用于绘制logistic模型
    theta:matrix[mx1]
    '''
    meshNum=50j
    xMesh,yMesh=numpy.mgrid[0:1:meshNum,0:1:meshNum]

    zMesh=numpy.zeros((50,50))
    for i in range(50):        
        for j in range(50):
            zMesh[i,j]=logistic_fun(xMesh[i,j],yMesh[i,j],theta)
            

    ax.plot_surface(xMesh,yMesh,zMesh,rstride=1,cstride=1,cmap=plt.cm.hsv,linestyles='dashdot',linewidth=0.5,alpha=0.3)
    plt.show()
	
def logistic_fun(x,y,theta):
    a=theta[0,0]
    b=theta[1,0]
    c=theta[2,0]
    return numpy.power(numpy.e,a*x+b*y+c)/(1+numpy.power(numpy.e,a*x+b*y+c))

def self_plot3d(x,y,z):
    '''
    自定义三维绘图模块
    x:matrix[nx1]
    y:matrix[nx1]
    z:matrix[nx1]
    x y z数据在使用前必须先进行转换  通过A1转换为一维array 
    否则plot函数将会使其为高维函数进行绘图
    '''
    
     
    x= x.A1
    y= y.A1
    z= z.A1

    x1=list()
    x2=list()
    y1=list()
    y2=list()
    z1=list()
    z2=list()
    
    for index in range(num):
        if z[index]==1:
            x1.append(x[index])
            y1.append(y[index])
            z1.append(z[index])
        else:
            x2.append(x[index])
            y2.append(y[index])
            z2.append(z[index])
   
    ax.scatter(x1,y1,z1,c='r',marker='o')
    ax.scatter(x2,y2,z2,c='b',marker='^')

    
if __name__=="__main__":
    newton_main()